function [varargout] = vtc_concat(targetfile, vtcs, varargin)
% vtc_concat  - concatenate VTCs
%
% FORMAT:       newvtc = vtc_concat(targetfile, vtclist [, options]);
%         or    newvtc = vtc_concat(targetfile, vtc1, vtc2, ... [, options]);
%
% Input fields:
%
%       targetfile  filename of VTC to write
%       vtclist     cell array with VTCs to concatenated
%       vtc1, ...   single VTCs to concatenated
%       options     1x1 struct with optional settings
%        .prt       filename of new protocol file, otherwise keep first
%        .trans     transformation to apply on each VTC on reading
%                   either of 'psc' or 'z'
%
% Output fields:
%
%       newvtc      object handle to newly written VTC

% Version:  v0.8a
% Build:    902503
% Date:     Oct-25 2009, 3:44 AM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   ~ischar(targetfile) || ...
    isempty(targetfile) || ...
   (~all(isBVQXfile(vtcs, 'vtc')) && ...
    ~iscell(vtcs)) || ...
    isempty(vtcs)
    error( ...
        'BVQXtools:BadArgument', ...
        'You must give a valid target filename and some VTCs.' ...
    );
end
if all(isBVQXfile(vtcs, 'vtc'))
    vtcc = cell(1, numel(vtcs));
    for vc = 1:numel(vtcs)
        vtcc{vc} = vtcs(vc);
    end
    vtcs = vtcc;
end
for n = 3:nargin
    if isBVQXfile(varargin{n - 2}, 'vtc')
        vtcs{end + 1} = varargin{n - 2};
    end
end
nv = numel(vtcs);
if nv == 1
    error( ...
        'BVQXtools:TooFewArguments', ...
        'You need at least two VTCs to concatenate.' ...
    );
end

% options ?
newprt = '';
transopt = 0;
if nargin > 2 && ...
    isstruct(varargin{end}) && ...
    numel(varargin{end}) == 1
    options = varargin{end};
else
    options = struct;
end
if isfield(options, 'prt') && ...
    ischar(options.prt)
    newprt = options.prt(:)';
end
if isfield(options, 'trans') && ...
    ischar(options.trans) && ...
   ~isempty(options.trans)
    switch (lower(options.trans(:)'))
        case {'psc', '%'}
            transopt = 1;
        case {'z'}
            transopt = 2;
    end
end

% check vtcs argument for filenames
loadvtcs = false(1, nv);
for vc = 1:nv
    if ischar(vtcs{vc}) && ...
        exist(vtcs{vc}(:)', 'file') == 2
        loadvtcs(vc) = true;
    elseif numel(vtcs{vc}) ~= 1 || ...
       ~isBVQXfile(vtcs{vc}, 'vtc')
        error( ...
            'BVQXtools:BadArgument', ...
            'Cell array must contain either VTC objects or filenames.' ...
        );
    end
end
if any(loadvtcs)
    vtctiosz = BVQXfile(0, 'transiosize', 'vtc');
    BVQXfile(0, 'transiosize', 'vtc', 65536);
    rlsame = BVQXfile(0, 'config', 'reloadsame');
    BVQXfile(0, 'config', 'reloadsame', true);
    for vc = find(loadvtcs)
        if ischar(vtcs{vc})
            try
                vtcfname = vtcs{vc}(:)';
                vtcs{vc} = BVQXfile(vtcfname);
            catch
                clearbvqxobjects(vtcs(loadvtcs));
                BVQXfile(0, 'config', 'reloadsame', rlsame);
                BVQXfile(0, 'transiosize', 'vtc', vtctiosz);
                error( ...
                    'BVQXtools:BVQXfileError', ...
                    'Error opening VTC ''%s''.', ...
                    vtcfname ...
                );
            end
        end
    end
    BVQXfile(0, 'config', 'reloadsame', rlsame);
    BVQXfile(0, 'transiosize', 'vtc', vtctiosz);
end

% check headers
vh = vtcs{1};
vs = size(vh.VTCData);
vvol = zeros(1, nv);
nvol = vs(1);
vvol(1) = nvol;
vtr = vh.TR;
if ~isempty(newprt)
    prts = cell(1, nv);
    try
        prts{1} = BVQXfile(prts{vc});
        if ~isBVQXfile(prts{1}, 'prt')
            error('BADORNOPRT');
        end
    catch
        if isBVQXfile(prts{1})
            prts{1}.ClearObject;
        end
        prts{1} = [];
        newprt = '';
    end
else
    prts = {};
end
for vc = 2:nv
    cs = size(vtcs{vc}.VTCData);
    if vtcs{vc}.Resolution ~= vh.Resolution || ...
        vtcs{vc}.TR ~= vtr || ...
        any(vtcs{vc}.BoundingBox.BBox(:) ~= vh.BoundingBox.BBox(:)) || ...
        any(cs(2:4) ~= vs(2:4))
        clearbvqxobjects(vtcs(loadvtcs));
        clearbvqxobjects(prts);
        error( ...
            'BVQXtools:BadArgument', ...
            'VTCs must match in spatial size/position and temporal settings.' ...
        );
    end
    nvol = nvol + cs(1);
    vvol(vc) = cs(1);
    if ~isempty(newprt)
        if vtcs{vc}.NrOfLinkedPRTs == 1
            prts{vc} = vtcs{vc}.NameOfLinkedPRT;
        elseif vtcs{vc}.NrOfLinkedPRTs > 1
            prts{vc} = vtcs{vc}.NameOfLinkedPRT{1};
        end
        if ~isempty(prts{vc})
            try
                prts{vc} = BVQXfile(prts{vc});
                if ~isBVQXfile(prts{vc}, 'prt')
                    error('BADORNOPRT');
                end
            catch
                if isBVQXfile(prts{vc})
                    prts{vc}.ClearObject;
                end
                prts{vc} = [];
            end
        end
    end
end
vvol = 1 + [0,cumsum(vvol)];
if any(isemptycell(prts))
    clearbvqxobjects(prts)
    prts = [];
end

% combine PRTs first
if ~isempty(prts)
    for vc = 2:nv
        prts{1}.Concatenate(prts{vc}, cumsum(vvol(1:(vc - 1))), vtr);
    end
    try
        prts{1}.SaveAs(newprt);
    catch
        warning( ...
            'BVQXtools:BVQXfileError', ...
            'Error saving new PRT: %s.', ...
            lasterror ...
        );
    end
    clearbvqxobjects(prts);
end

% construct transioobject
tobjname = tempname;
tobj = transio(tobjname, 'ieee-le', 'uint16', 0, [nvol, vs(2:4)], 1);

% copy first VTC object, and set new NrOfVolumes, VTCData
newvtc = bless(vtcs{1}.CopyObject, 1);
newvtc.NrOfVolumes = nvol;
newvtc.VTCData = tobj;

% create temporary object in mem
mo = uint16(0);
mo(nvol, vs(2), vs(3)) = 0;

% loop over slices
for sc = 1:vs(4)
    for vc = 1:nv
        if transopt == 2
            mo(vvol(vc):(vvol(vc + 1)-1), :, :) = ...
                16000 + 4000 * ztrans(double(vtcs{vc}.VTCData(:, :, :, sc)));
        elseif transopt == 1
            mo(vvol(vc):(vvol(vc + 1)-1), :, :) = ...
                100 * psctrans(double(vtcs{vc}.VTCData(:, :, :, sc)));
        else
            mo(vvol(vc):(vvol(vc + 1)-1), :, :) = ...
                vtcs{vc}.VTCData(:, :, :, sc);
        end
    end
    newvtc.VTCData(:, :, :, sc) = mo;
end

% clear objects
clearbvqxobjects(vtcs(loadvtcs));

% save VTC under new name
try
    if ~isempty(newprt)
        newvtc.NameOfLinkedPRT = newprt;
    end
    newvtc.SaveAs(targetfile);
catch
    warning( ...
        'BVQXtools:BadArgument', ...
        'Bad target filename given. Please remove ''%s'' manually.', ...
        tobjname ...
    );
    newvtc.ClearObject;
    return;
end

% delete old temp object
delete(tobjname);

% if no argout, clear object
if nargout < 1
    newvtc.ClearObject;
else
    varargout = cell(1, nargout);
    varargout{1} = newvtc;
end
